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Fourier transform scanning tunneling spectroscopy (FT-STS) measures the scattering of conduc¬ 
tion electrons from impurities and defects, giving information about the electronic structure of both 
the host material and adsorbed impurities. We interpret such FT-STS measurements in terms of 
the quasiparticle interference (QPI), here investigating in detail the QPI due to magnetic impurities 
adsorbed on a range of representative non-magnetic host surfaces, and contrasting with the case of 
a simple scalar impurity or point defect. We demonstrate how the electronic correlations present for 
magnetic impurities markedly affect the QPI, showing e.g. a large intensity enhancement due to the 
Kondo effect, and universality at low temperatures/scanning-energies. The commonly-used joint 
density of states (JDOS) interpretation of FT-STS measurements is also considered, and shown to 
be insufficient in many cases, including that of magnetic impurities. 
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I. INTRODUCTION 

Magnetic atoms embedded on non-magnetic surfaces 
provide realizations of quantum impurity models^ which 
are amenable to detailed experimental study and manip¬ 
ulation with a scanning tunneling microscope (STM).^“^ 
Such systems are promising candidates as a basis for 
nanoscale computational, memory storage and spintronic 
devices.They are also of fundamental interest in their 
own right, due to the subtle interplay of strongly corre¬ 
lated local spin and orbital degrees of freedom coupled 
to a conduction electron bath. 

In scanning tunneling spectroscopy (STS) experi¬ 
ments, the differential conductance between tip and 
surface is related to the local density of electronic 
states (LDOS) at a particular scanning energy and 
temperature;^’^ by rastering the STM tip across the sur¬ 
face (Fig. 1) a spatial map of the LDOS in the vicinity of 
features such as adsorbed impurities may be generated.^ 
Such impurities break translational symmetry at the sur¬ 
face, causing scattering of conduction electrons and mod¬ 
ulations in the LDOS that depend strongly on both the 
electronic structure of the underlying sample (‘host’) and 
the properties and distribution of impurities. Fourier 
transform STS (FT-STS) - in which such modulations 
are analysed in reciprocal space and interpreted in terms 
of the quasiparticle interference (QPI) between the diag¬ 
onal states of the clean host^^^^ - thus provides a wealth 
of information on the nature of the impurity-host and 
(host-mediated) inter-impurity correlations. 

The aim of this work is to calculate the QPI due to 
magnetic adatoms on metallic surfaces, drawing compar¬ 
ison with the case of non-magnetic, ‘scalar’ impurities 
(s-wave potential scatterers in the weak-scattering Born 
limit^^). We also link quantitatively FT-STS measure¬ 
ments to the calculated QPI for generic systems. By 
way of contrast, we then examine the widely-used phe¬ 
nomenological joint density of states (JDOS) interpreta¬ 



FIG. I. Schematic FT-STS setup: a spatial map of surface 
LDOS modulations due to electron scattering from a magnetic 
adatom impurity is extracted from differential conductance 
measurements as an STM tip is rastered over the surface. 


tion for FT-STS measurements. 

In addition to FT-STS, magnetic impurities can them¬ 
selves be spectroscopically probed when the STM tip is 
positioned directly over an impurity. A narrow Kondo 
resonance is commonly observed around the Fermi en¬ 
ergy in these local STS measurements,^’^ arising from 
the formation of a many-body Kondo singlet state where 
the impurity local moment is screened dynamically by 
surrounding host conduction electrons. For STS stud¬ 
ies of Go adatoms on noble metal surfaces, typical 
Kondo temperatures of Tk ^ 50 — lOOK are extracted 
from the half-width at half-maximum of the spectral 
Kondo resonance.The single-orbital Anderson impu¬ 
rity model^^ has been successfully used to rationalize 
such local STS measurements.It is adopted here 
as a qualitatively accurate description of a generic, low- 
spin magnetic adatom; and is treated using the numeri¬ 
cal renormalization group (NRG).^^ In other cases, e.g. 
high-spin Mn or Fe adatoms, generalized impurity mod- 
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els must be employed to capture the full orbital structure 
of the adatom aud the host material^^’^^ (such calcula- 
tious are quite specific to the particular system uuder 
cousideratiou,^^ aud are uot explicitly cousidered). 

We begiu by cousideriug a geueral formulatiou 
(Secs. II,III) for auy uumber or type of impurities em¬ 
bedded ou the surface of tight-biudiug hosts. The siu- 
gle maguetic impurity case, which is the maiu focus of 
this paper, is simply a particular example (we subse- 
queutly apply the formalism to multi-impurity systems 
iu Ref. 23). We highlight iu particular (Sec. Ill) the dif- 
fereuces betweeu the 3d cubic lattice with a (100) surface, 
the 2d square lattice, aud the houeycomb lattice. These 
simple but represeutative hosts reproduce a rauge of pos¬ 
sible material realizatious: the LDOS for each features 
distiuct behavior close to the Fermi level, with flat, di- 
vergeut aud vauishiug LDOS for the three lattices respec¬ 
tively. Each host lattice gives rise to qualitatively differ- 
eut impurity physics, with e.g. the siugle-impurity Koudo 
temperature siguificautly euhauced (suppressed) by iu- 
creased (depleted) deusity of states arouud the Fermi 

level.24-27 

Siugle impurity systems are cousidered explicitly iu 
Sec. IV ff. The QPI due to siugle maguetic aud scalar im¬ 
purities for each host is cousidered iu Sec. IV, drawiug at- 
teutiou to the qualitative differeuces iu QPI arisiug from 
the differeut types of scatteriug ceuter aud host. The 
rich dynamical properties of the QPI are studied iu de¬ 
tail for systems coutaiuiug a maguetic impurity iu Sec. V; 
we emphasize that this stroug euergy-depeudeuce cauuot 
be reproduced iu systems coutaiuiug ouly scalar impuri¬ 
ties or structural defects. ludeed, stroug electrou corre- 
latious iu systems with maguetic impurities are showu to 
produce uuique QPI siguatures: iu metallic systems, the 
QPI exhibits uuiversality iu terms of rescaled scauuiug 
euergy aud temperature due to the Koudo effect, while 
uou-trivial local momeut physics iu observed iu the 2d 
houeycomb case.^^’^^’^^“^^ 

Sec. VI examiues the iuterpretatiou of experimeutally- 
measurable FT-STS. By simulatiug the experimeutal 
protocol, we iuvestigate the possible deviatiou of FT-STS 
measuremeuts from calculated QPI due to the fiuite size 
of the LDOS plaquette measuremeuts iu real-space. We 
couclude by critically examiuiug the relatiouship betweeu 
the JDOS aud the QPI, showiug that the two quautities 
are uot directly related iu auy but the simplest case (a 
siugle scalar impurity embedded ou a ceutrosymmetric 
host), where the JDOS aud QPI are Hilbert coujugates. 


II. HOST SYSTEMS IMPURITY PROBLEM 
A. Model 

We cousider a host material with impurities deposited 
ou the surface, which scatter the quasiparticles of the 
cleau host. luitially we cousider the geueral problem of N 
maguetic surface adatoms, formulatiug the QPI geuerally 


for auy type, uumber aud distributiou of adatoms. 

The cleau host is takeu to be uou-iuteractiug aud of 
tight-biudiug form, giveu iu its real-space basis by 

"Hhost = eo F “ t F 

1 , 0 - {ij),o- 

where creates au electrou of spiu a =t/i iu the Wau- 
uier orbital localized at site r^; aud (ij) deuotes a sum 
over uearest-ueighbor sites, coupled by a tuuuel matrix 
elemeut t. Here we cousider the half-filled host, 6 q = 0. 

Specifically, we focus ou the simple 2d square lattice, 
the 3d cubic lattice with a (100) surface, aud the 2d hou¬ 
eycomb lattice. The 3d cubic lattice iu particular is rep¬ 
reseutative of a wide class of regular metallic systems, 
with a coustaut (fiuite) electrouic deusity of states at 
low euergies. By coutrast, the 2d square lattice, relevaut 
to certaiu layered materials, features a vau Hove siugu- 
larity with a logarithmically divergiug deusity of states 
arouud the Fermi level.The houeycomb lattice, de- 
scribiug grapheue withiu the simplest uou-iuteractiug 
tight-biudiug approximatiou, is uotable because it is a 
bipartite lattice with a vauishiug (pseudogapped) deusity 
of states at low euergies.Although simplified, these 
host systems exemplify a uumber of distiuctive features 
relevaut to real materials, aud each iuduces qualitatively 
differeut impurity physics (Sec. IV). 

The full model, iucludiug N impurities, is giveu by 

N 

^ — ^host T E (2) 

a=l 

Iu the simplest case, the local poteutial at site is 
modified by impurity a, breakiug trauslatioual iuvariauce 
aud causiug additioual electrouic scatteriug. These static 
defects, referred to as ‘poteutial scatteriug’ (or ‘scalar’) 
impurities, are described by 

(3) 

a 

This simple model, while ofteu appropriate to describe 
poiut defects iu materials, does uot faithfully capture 
the physics of mauy adsorbed impurity adatoms — iu 
particular maguetic impurities, which are dynamic ob¬ 
jects with iuterual degrees of freedom aud stroug local 
Coulomb iuteractious.^ These have a siguificaut qualita¬ 
tive effect ou the scatteriug properties, aud must there¬ 
fore be takeu iuto accouut. Iu this paper we focus primar¬ 
ily ou adsorbed maguetic impurities, described iu terms 
of correlated quautum levels 

qjAIM _ f n 4- d d^ d 

" (4) 

+y^(dLcr^, + H.c.) , 

a 

where d]^^ creates a spiu-a electrou ou impurity a, which 
is coupled to host site by tuuueliug matrix elemeut 
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V. For simplicity we consider explicitly the particle-hole 
symmetric case €d = —1//2. We emphasize that true 
magnetic impurities preserve SU(2) spin symmetry and 
time-reversal symmetry, unlike a static polarized local 
moment, which simply acts like a local magnetic field. 


B. Impurity dynamics 

Single-particle dynamics of correlated impurities em¬ 
bedded in a non-interacting host are described generically 
by the Green function matrix [Gd{^)]a,/3 = {u) = 

with {{A]B))^ the Fourier transform of 

the retarded correlator —i0{t){{A{t),B{0)}). The ma¬ 
trix Dyson equation for the impurity Green functions is 

[Grf(c^)]-' = [grfH]-' - SH . (5) 

The non-interacting (but host-coupled) impurity propa¬ 
gators are given by 

[gd(w)]“^ = (w + z 0 +- ed)I - r(w) , (6) 

in terms of the hybridization matrix T{uj) with elements 
[r(w)]„,;3 = V'^G°{ra,rp,uj); where G°(ra,r/3,w) = 
is the propagator between sites and 
of the clean host (without impurities), and and are 
the host sites to which impurities a and [3 are coupled. 
The self-energy matrix contains all the nontrivial 

information due to electronic interactions, which give rise 
to the Kondo effect, RKKY interaction, etc. 

In the present work we employ NRG to solve the un¬ 
derlying quantum impurity problem. In the spirit of 
Ref. 34, equations of motion can be used to obtain an 
expression for the self-energy matrix, 

S(a;) = [GdH]-iFdH , (7) 

where [Fd(w)]„,/3 = U Both Gd(w) 
and Fd{oo) are calculated directly in NRG using the 
full density matrix approach^^’^^ within the complete 
Anders-Schiller basis. 

In the case of a single impurity, the Dyson equation 
(Eq. 5) reduces to Gy{uj) = (c(;+iO+— —r(c(;) — En)”^, 

with the hybridization r(co’) = 'F^G^(ri, ri, cj) related 
simply to the clean host LDOS (itself independent of po¬ 
sition due to translational invariance). In NRG,^^ a dis¬ 
cretized version of the conduction electron Hamiltonian 
is formulated, and mapped onto a semi-infinite Id chain 
with the impurity located at one end. Discretizing on a 
logarithmic energy scale leads to the ‘Wilson chain’ rep¬ 
resentation in which hopping matrix elements decrease 
exponentially down the chain. The RG scheme involves 
iterative diagonalization, starting at the impurity and 
working down the Wilson chain, discarding high-energy 
states at each step.^^ 


C. Host dynamics 

The full dynamics of the host, in the presence of 
impurities, is embodied in the lattice Green functions 
G(r^,rj,co’) connecting arbitrary host sites and Vj. 
They can be related exactly to the above impurity Green 
functions by equations of motion: 

G(ri,rj,w)-G°(rj,rj,w) = 

'^G°{ri,rc„ui)Tap{ui)G°{rp,rj,ui) , (^) 

a,(3 

where the sum runs over impurities a and and Tagioo) 
is the real-space t-matrix. For magnetic impurities de¬ 
scribed by in Eq. 4, the t-matrix takes the form 

= V^Gfiuj) , ( 9 ) 

requiring as such a knowledge of the full impurity Green 
functions from Eq. 5. By contrast, the t-matrix for po¬ 
tential scattering impurities, described by in Eq. 3, 
can be obtained simply in closed form, 

Trpi^)=v[l-vG^iu^)]-^^, (10) 

where the elements [G^(co’)]c,'/ 3 ' = G^(rc,',r/ 3 /,cj) are free 
host Green functions. 

In the diagonal quasiparticle basis of the clean host, 

^host = ^^^k^kcr^kcr 
k,cr 

with Ck the dispersion (k labels the Bloch state momen¬ 
tum). The t-matrix equation (Eq. 8 ) can also be trans¬ 
formed into the momentum-space basis, 

AG(k, k', uj) = G^(k, uj)T{k, k', a;)G^(k', u) , ( 12 ) 

where AG(k, k',C(;) = G(k, k',^^) — G^(k, C(;)^k,k' and 
G^(k, cj) = (cj + i0+ — ek)“^. All quasiparticle scattering 
induced by the impurities is now contained in 

T(k,k',w) = ^ X To^pito) , (13) 

where is the volume of the first Brillouin zone. 

III. QPI 

At sufficiently low temperature and bias, the differen¬ 
tial conductance (i/(r, V)jdV between STM tip and sur¬ 
face for a particular bias V is proportional to the LDOS, 
p(r,co’ = eE), of the sample at position r. STS thus 
measures the energy-resolved electronic structure of the 
sample at a particular point in real space. 

The experimental quantity of interest is the QPI, 

P(q>w)= E e"“i"'V(ri,w) , (14) 

ie(LxL) 
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obtained at given scanning energy uj. It is therefore the 
Fourier transform of a real-space LDOS map p(ri,co’), it¬ 
self measured by STS over an L x L sample region. Since 
the STS experiment probes the surface, the LDOS map is 
two-dimensional. In 2d systems such as graphene, or ef¬ 
fective 2d layered systems such as the cuprates, STS thus 
probes directly the underlying lattice. In 3d systems, by 
contrast, the lattice probed by STS corresponds to the 
crystallographic surface lattice. 

The FT-STS technique has been employed to inves¬ 
tigate the electronic structure of materials and to map 
their Fermi surface contours,providing complemen¬ 
tary information to techniques such as ARPES.^^’^^ FT- 
STS has proven of particular value in the study of layered 
materials such as cuprates and pnictides,^^™^^ as well as 
topological insulators,graphene^^ and heavy fermion 
materials.^^ In many cases (e.g. for weak, disorder- 
induced scattering in cuprates^^), the local defects giving 
rise to the QPI may be approximated as scalar impuri¬ 
ties in the Born limitalthough even for scalar impuri¬ 
ties this simplification is known to be insufficient in some 
cases.For e.g. transition metal adatoms, however, the 
full interacting impurity model must be considered. 

The QPI at q = 0 is often omitted in experimental re¬ 
sults, because it corresponds to the total density of states 
sampled, and is extensive in L^. The desired impurity 
contribution to the QPI is then obtained by subtracting 
the result for the clean system without impurities, 

Ap{q,Lo)= (15) 


lattice with any number or type of impurities. Although 
the LDOS is sampled over a finite region, we stress that 
the host system is in the thermodynamic limit. 

The accurate calculation of lattice Green func¬ 
tions G^{Ti,Tj,(jo) is itself a subtle and well-studied 
problem.Exact diagonalization of finite-sized lat¬ 
tices or discrete Eourier transforms yield poor approxi¬ 
mations to Green functions of the desired (semi-)infinite 
systems, especially at low scanning energies or near van 
Hove singularities. Recursion methods^^’^^ can be used 
for periodic systems where the exact dispersion ek is 
known — but such techniques are numerically unstable 
for large site separations |r^ —r^j, and converged solutions 
can be computationally demanding.Although ‘bond 
cutting’^^ and ‘continued fraction’^^ variants have been 
developed in special cases, recursion methods typically 
cannot be used for systems with a surface that breaks 
translational invariance, limiting applicability within the 
QPI context. In consequence, calculation of all N x 
non-local Green functions required for a system of N im¬ 
purities in an L X L region is often the major challenge. 

In this paper we have developed instead a novel tech¬ 
nique for fast and accurate numerical calculation of free 
Green functions on hypercubic-type lattices. The method 
is detailed in the appendix, and involves successive con¬ 
volutions of simpler Id Green functions which are known 
exactly in closed form. The convolution itself can more¬ 
over be performed efficiently using fast Eourier transform. 

B. Scattering state formulation 


where Ap = p — (with p^ for the clean host). 
Since p^(q, cj) (x ^^(q^^) scales with the num¬ 

ber of impurities. The normalized QPI power spectrum, 

is intensive, indepen¬ 
dent of the number of impurities or sample region size. 

We now give a general formulation for calculating the 
QPI due to scattering from single or multiple impurities, 
which can either be simple static potential scattering de¬ 
fects or magnetic (dynamic, interacting) impurities. 


A. Real-space formulation 

Eollowing the experimental protocol, the QPI can be 
calculated by discrete Eourier transform of the LDOS 
within an 1/ X 1/ region of the host surface, using Eq. 15. 
The LDOS at site in the presence of impurities is 
related to the local host Green function, p{Ti,uj) = 
— ;^Im G{Ti,Ti,(jo), such that from Eq. 8, 

Ap{ri,ui) = -tlm y'G°(ri,ra,w)Ta/3{w)G°(r^,ri,w) 

TT ^ 

a,(3 

(16) 

in terms of the full scattering t-matrix and free non-local 
host Green functions. The virtue of Eq. 16 is that it 
is entirely general, and can in principle be used for any 


The ‘true’ QPI containing all scattering information is 
obtained by taking the thermodynamic limit of plaque- 
tte size, L ^ oo. Translational symmetry implies a basis 
of states with well defined momentum parallel to the sur¬ 
face, i.e. over the first surface Brillouin zone (ISBZ). The 
local Green function for a surface lattice site is expressed 
in terms of this basis by 2d Eourier transformation. 


G{ri,ri,u)) = JJ 


d^kiid^kJi 

^BZ 




xG(k||,kj|,cj) , 


ISBZ 


(17) 

with VIbz the volume of the ISBZ. Writing Ap{ri,uj) = 
— ;^Im AG{ri,ri,uj), Eq. 15 takes the form. 


Ap(q,c^)=-i 

TT ^ 
i 

d^kiid^kn 


X Im 


.-irr(k;|-k||) 


n 


ISBZ 


BZ 


AG(k||,kj|,a;) , 
(18) 


where AG(k||,k',w) = G(k||,k',w) - G°(k||,k',w). 
Eq. 18 can itself be recast as 

Ap(q, co) = -f^ [Q(q, uj) - Q(-q, w)*] , (19) 
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where 

Q(q>‘^)= j AG(k||,k|| - q,w) ( 20 a) 

ISBZ 

= '^Tapioj) X Aa/ 3 (q,w) . (20b) 

Ck:,/3 

As highlighted by Eq. 20 b, the QPI factorises into 
a momentum-independent scattering amplitude Tc^(co’) 
(Eqs. 9 or 10), and a host response function Ac^^(q, cj) 
which depends only on the host lattice and the spatial 
location of impurities, but not the type of impurity (and 
thus details of the scattering). The explicit form of this 
host function must be determined separately for each lat¬ 
tice, as considered below. 


1. 2d square lattice 

Consider first the 2 d square lattice, where the QPI 
calculation is simplest. As the system is itself two- 
dimensional, the surface-momentum basis is simply the 
diagonal representation, k = ky. The QPI thus follows 
from Eqs. 20a, 12,13, and is indeed of form Eq. 20b with 


plane in isolation - but surface states labelled by ky re¬ 
main coupled to the bulk (and thus to each other). In 
general, AG(ky,ky — q,cj) does not therefore take the 
form of Eq. 12, but rather 


AG(k||,kj|,w) = 

G°(k||, k||, w)T(kj|, kj",w)G°(k||', k||, w) , 

ISBZ 

(23) 

where G^(ky, kj|, cj) is a complex through-bulk propaga¬ 
tor. Calculation of the QPI thus in general requires the 
full integrals over intermediate scattering pathways. 

In the case of hypercubic-type lattices, however, a sig¬ 
nificant simplification arises, because ky is still a good 
quantum number; so a surface momentum t-matrix equa¬ 
tion with the same structure as Eq. 12 still applies, albeit 
with modified host surface Green functions. The result 
for the 3d cubic lattice with a (100) surface is simply, 

AG(k||,kj|,c^) =GLf(k||,w)T(k||,kj|,c^)GLf(kj|,c^) , 

(24) 

where the t-matrix is still given by Eq. 13, but 


f f d^k''d^k'' 

J J ^BZ 


Ac«/3(q,w)= [ ^G0(k,a;)G°(k-q,a;)e‘[‘^--“-(‘^-‘i)--'’l 

J ^tBZ 

ISBZ 

( 21 ) 

The free momentum-space Green functions are them¬ 
selves given by^^ G^(k,C(;) = {uj i 0 + — ek)“^, with 2 d 
square lattice dispersion (and lattice constant ao) 

ek = — 2 t[cos(aoka,) + cos(aok^)] . ( 22 ) 

The half-bandwidth is then D = At in terms of the lattice 
hopping matrix element t appearing in Eq. 1. 

Aa/ 3 (q, ^) can be computed efficiently by using the con¬ 
volution theorem to do the Brillouin zone integration: 

-Fk[A„^(q,w)] = Jk[G0(k,:^)e‘*'--“] X 

where Jdc denotes the 2d fast Eourier transform. 


2. 3d cubic lattice with (100) surface 


GLf(k||,a;)=/ 



where 


tf{ui) =io- 


sgn(A)\/A2 — 1 

iVl — 


|A| > 1 
|A| < 1 


(25) 


with ek|| the 2d square lattice dispersion, Eq. 22. This 
broadens the pole in G^(k, cj) arising for the pure 2 d 
system, to an ellipse of width 2t centered on uj = eky 
in Gg^^£(ky, cj) for the 3d system. It follows that the 
structure of the host function Aap{ci,uj) is the same as 
in the 2d square case, Eq. 21, with Gg^^f(k||, cj) in place 
of G^(k,co’). In this paper we calculate the true 3d cubic 
QPI via Eqs. 19, 20, 21 , using the exact expression for 
the bulk-coupled surface Green functions, Eq. 25. 

Eor more complex systems where such a prescription 
is not available, the bulk dephasing of pure 2 d surface 
states could be approximated by using G^(ky,co’) = (cj + 

~ Ck||)“^, with finite > 0. Green function poles are 
thereby lifetime broadened by lorentzians of width 77 . 


3d host lattices are more subtle, due to the surface- 
sensitive STM measurement. As only the surface LDOS 
is probed, the QPI amounts to a partial trace over the 
t-matrix equation Eq. 12, in contrast to the full trace for 
the 2d square lattice (Eq. 21 ). The QPI must thus be 
evaluated in a basis which preserves the layer index, the 
surface momentum basis of Eq. 17 (rather than the diag¬ 
onal basis of Eq. 12 ). AG(ky,ky — q,cj) in Eq. 20 thus 
involves propagators between states with surface momen¬ 
tum ky and ky — q. The 2 d transform of G(r^, r^, cj), 
Eq. 17, leads to a diagonal representation in each 2d 


3. Honeycomb lattice 

The honeycomb lattice is complicated by the bipartite 
nature of the lattice, which is viewed as two interlocking 
triangular sublattices. We define 

t(k) = [1 + (26) 

in terms of the triangular sublattice vectors ai and a 2 
and the inter-sublattice vector r = — rf (where rj is 
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a site i on the ^ = AjB sublattice). The A/5-sublattice 
structure gives rise to distinct +/— bands in momentum 
space.The honeycomb lattice dispersion for these +/— 
bands, and the complex phase, follow as 

d = ±l^(k)l (27a) 

0(k) = arg [t(k)] . (27b) 

Real-space operators are expressed in a diagonal basis by. 


C 7 — 


1 




IBZ 




1/2 

BZ 


gir/’k ^is^(/)(k )/2 


^-,k + '^7^+,k 


where 5 ^ = ±1 for the A or B sublattice. 


(28) 


Generalizing Eqs. 15 and 16 to take account of this 
sublattice structure gives 


Ap(q,w) = E ^Pi,71,72 5 


(29) 


7^71,72 


where 




aeii (3^72 

(30) 


with the real-space sum over i spanning sites rj on sub¬ 
lattice 7 . Impurity a{P) is taken to be on sublattice 
71 ( 72 )- Thus, is contribution to the 

full QPI from sites on the 7 sublattice due to impurity- 
induced scattering between 71 and 72 sublattices. 

Using Eq. 28 in the definition ,uj) = 

((c^^; ))^, takes the same form as 

Eq. 18,' 

TT 

X Im JJ , 

IBZ 

(31) 

but now with 

^<^7,7l,72(k,k',w) = E E AbZ 

aeii P^72 

X 0i[(s7-«7l)^(k)-(s7-«72)'^(^')]/2 

X (s^GO (k) + s^.G^k)) (s^GKk') + s^,G7k')) , 

(32) 


in terms of the +/— band free Green functions^^ 

Gl(k) ^ ((c±,k;4,k))° = (^ + i0+ - 4)“'- The QPI 
contribution then follows as. 


^^ 7 , 71,72 (^5 [^ 7 , 71,72 (^5 Q 7 , 7 i, 72 ( Q 5 ] 

(33) 


where 

e7,7i.72(q>‘^) = E E Bp{oj) X A7'’^"(q,w).(34) 

aeii /3e72 

Eqs. 33 and 34 are thus analogues of Eqs. 19 and 20, with 

[ ^ ^ p-i(rT-k-rp-(k-q)) i[(g^-g^ 70 (k)-(g^-g^^) 0 (k-q )]/2 

J ^^BZ 

IBZ 

X ( 5 ^G^(k) “h 5 ^^G^(k)) ( 5 ^G^(k — q) + ~ q)) 

(35) 

We stress that interband scattering and the momentum- 
dependent phase factors appearing in Eq. 35 are impor¬ 
tant, and affect the full QPI qualitatively. 


IV. SINGLE IMPURITY QPI: 
MAGNETIC SCALAR IMPURITIES 

The generalized problem involving N magnetic im¬ 
purities, spatially separated and coupled to conduction 
electrons of the host lattice, is naturally highly rich and 
complex (we consider aspects of it in subsequent work^^). 
Erom here on in this paper, we focus on a single magnetic 
impurity - in terms of which experimental QPI patterns 
are in fact typically interpreted. In this case the QPI is 
given by Eq. 20b (with a = (3 = dropped hereafter), 

Q(q, Lo) = T(w)A(q, w) , (36) 

with T{uj) the single-impurity t-matrix and explicit forms 
for the host functions A(q, cj) given in Secs. IVB-IVD. 

The impurity itself is often taken to be a static po¬ 
tential defect in the weak-scattering Born limit.The 
t-matrix is then pure real and energy-independent, 

T^^{UJ) 73 V , (37) 

with V the potential scattering strength (see Eq. 3). 
Eq. 37 is the leading-order approximation to the exact 
Eq. 10, holding provided \vG^{yo^^Y(^^uj)\ ^ I. We add 
however that this approximation is not valid in the vicin¬ 
ity of divergences in the host density of states (arising e.g. 
at cj = 0 in the 2d square lattice). In the special case of a 
single impurity on a centrosymmetric surface, ( 3 (q, cj) = 
(5(—q,cj), so (Eq. 19) Ap(q,cj) = —(5(q,cj). Eor a 
scalar impurity in the Born limit, 

Ap(q,c^) --A"(q,c^) (38) 

TT 

where A(q, uj) = A'(q, cj) +iA"(q, cj). The QPI scanning- 
energy dependence is thus due solely to that of A"(q, cj). 

Eor magnetic impurities by contrast, electron correla¬ 
tions give rise to nontrivial dynamics. Erom Eq. 9, 

r'"“9(w) = V^Gdiuj) 


(39) 
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2D Square 3D Cubic (100) 2D Honeycomb 



FIG. 2. Dynamics of a single magnetic impurity on various lattices at T = 0, plotted as —Im tGd{(^) (solid lines) and Re tGd{cj^j) 
(dashed) \s cu/t, calculated via NRG. Insets show corresponding host density of states tpo(<^)- Impurity parameters: U = 1.95t 
and V = 0.555t for 2d square and 3d cubic lattices; U = 0.704t and V = 0.493t for honeycomb lattice. 


in terms of the impurity Green function Gd{oo); the QPI 
follows from Eqs. 19,36. For a centrosymmetric surface, 

A/9(q, w) Gd(w)A"(q,a;) 

+ Im G'd(w)A'(q,w)] , 

with contributions from both real and imaginary parts of 
A(q, cj), and weights that depend on the impurity Green 
function at energy uj. As discussed below, the Kondo 
effect produces a scattering enhancement at low temper¬ 
atures and scanning energies, causing a crossover in the 
QPI from being dominated by A"(q, cj) at high energies 
(similar to that of a scalar impurity) to being dominated 
by A'(q, cj) at low energies. 


A. Effect of host on impurity dynamics 


For a magnetic impurity, the QPI depends on both the 
host function A(q, cj) and the impurity Green function 
Gd{^) - which itself depends on the host. Specifically, 
the impurity problem is controlled by the hybridization 
function F(co’), related to the clean host density of states, 
Po(<^), by —Im r(co’) = 7rV‘^po{uj). The Kondo physics 
is sensitive to the behavior of po{uj) near the Fermi level 
(cj = 0);^^“^^’^^ and to leading order. 


Po{^ 




log( 16 t/|u;|) 

27r‘^t 

J_ _ 

6\/27r^t3 

1^1 


2d square 

3d cubic, (100) surface 
honeycomb . 


(41) 

These lattices exemplify three paradigms, with densities 
of states that are diverging, flat, or pseudogapped at low- 
energy. po{uj) vs uj/t is shown in the insets to Fig. 2. 

The density of states for metallic systems is typically 
flat at low energies. This gives rise to an exponentially- 
small Kondo scale^ Tx/t ^ exp[—7rI7/8K^po(0)], and 


low-energy Fermi liquid behavior^ 

ImF(cj) X lmGd{oj) 1 - a^i^ujjTKf + ... , (42) 

with Fermi level spectral pinning, —7rV‘^po{0)lmGd{0) = 
1 (as ImF(O) = —7rV^po{0)). This is shown for a mag¬ 
netic impurity on the (100) surface of a 3d cubic lattice in 
Fig. 2 (center panel), where the imaginary and real parts 
of tGd{oj) are plotted vs uj/t. We have chosen represen¬ 
tative impurity parameters U = 1.95t and V = 0.56t, 
yielding Tjc ~ 5 x 10“^t [defined here as the half-width 
at half-maximum of the Kondo resonance^^]. With the 
host bandwidth 12t = lleV (such that U = 1.79eV and 
V = 0.51eV) we obtain Tk ~ 57K, consistent with es¬ 
tablished results^^47,56 atoms on a Cu surface. As 

seen from Fig. 2 (center), the Kondo effect results in a 
large imaginary part —tim Gdioj) ^ 6t^/(7rK^) 6.2 for 

low energies |ci;| ^ T^. 

For the 2d square lattice, the low-energy divergence in 
the host density of states results in an enhanced Kondo 
temperature.^^ In any Kondo phase, the pinning condi¬ 
tion from Eq. 42 still holds,implying that, 

-lmGd(w) ^[ln(16t/|w|)]“h (43) 

which decays logarithmically at low energies. As con¬ 
firmed in the left panel of Eig. 2, the impurity spectrum 
—Im tGdioj) therefore shows a maximum at \uj\ ^ Tk- 
With the same parameters as the 3d cubic system, we 
now obtain a much higher Tk ~ 584K.^^ 

Einally, in the pseudogapped honeycomb lattice the 
Kondo effect is suppressed due to the depleted density of 
states near the Eermi level, and the local moment phase is 
stable for any JJjV^ at particle-hole symmetry.The 
impurity spectrum then takes the low-energy form 

-IinGrf(w) a'\Lu\ (44) 

(with a' a constant). This decay of Gd{oj) is demon¬ 
strated in the right panel of Eig. 2. With t = 2.84eV, the 














































FIG. 3. QPI for a single impurity on the 2d square lattice at 
a scanning energy uj = 0.055t ^ Tk- Upper panels compare 
the QPI maps |Ap(q)| for scalar (left) and magnetic (right) 
impurities; lower panel shows a Brillouin zone cut along the 
path r ^ X ^ M ^ r, plotted as Ap(q)/Aptot. Symmetry 
points defined as qr = 0 , qx = Ai, qM = A 1 +A 2 in terms of 
reciprocal lattice vectors Ai = 27r/ao(l, 0), A 2 = 27r/ao(0,1). 
Magnetic impurity parameters as in Fig. 2 such that uj — 
5 X = Tk] and v — 0.5t for the scalar impurity. 

honeycomb lattice models the tt/tt* bands of graphene; 
we use U = 0.704t and V = 0.493t as realistic impurity 
parameters obtained from ab-initio calculations for Co 
atoms on graphene. 

We now turn to the QPI for these three lattices. 


B. 2d square lattice 

The QPI is obtained from Eq. 38 or 40, with A(q, uj) for 
a single impurity given from Eq. 21 (with Tq, = = 0), 

/ rj^ Ic 

-— G°(k, uj)G°{k - q, w) . (45) 

ISBZ 

Eig. 3 shows the absolute value of the QPI |Ap(q)| 
as a colour map in q-space (upper panels), compar¬ 
ing the scalar impurity (left) with the magnetic impu¬ 
rity (right), at a fixed scanning energy uj Tk, using 
the same parameters as Eig. 2. The lower panel shows 
a cut across the Brillouin zone of Ap(q)/Aptot 5 where 
Aptot = Aptot(<^) = 

tering amplitude at energy uj. Eor a single impurity, the 
topology in q-space, on which we now focus, is completely 
determined by the host function A(q, uj) (see Eqs. 38, 40). 

Eor scalar impurities in the Born limit, only the imag¬ 
inary part, A"(q,cj), plays a role (Eq. 38). Eor given 
scanning energy cj, its structure gives rise to singular 
lines at q = q*(ci;) associated with the van Hove point 


FIG. 4. QPI for an impurity on the (100) surface of a 3d cubic 
lattice. As Fig. 3 but at a; = 5 x 10“^t = Tk- 

of the 2d square lattice. Around the F symmetry point, 
these lines are rectangular hyperbola, 

±ar(a;)2 : -n < < n , (46) 

with the dispersive properties controlled by ar(ci;); while 
around the M point 

(27r-aM(w))^ : - ir < q*^y < ir . 

(47) 

In q-space, A"(q, cj) is found to diverge as 

A"(q,a;) |q-qT'^' (48) 

when approaching a point q* from the F or M points. The 
region enclosed by these divergences is therefore charac¬ 
terized by high QPI scattering intensity — see Fig. 3 for 
the scalar impurity. A"(q, cj) does not however diverge 
on approaching from X, and remains comparatively small 
in its vicinity. In fact, A"(q, uj) is odd in uj due to the ex¬ 
act symmetry A(q, cj) = A(q, —u;)*. The scalar impurity 
QPI Ap(q, cj) ^ UJ thus vanishes at low energies away 
from the lines of divergence. 

The situation is rather different for the magnetic im¬ 
purity because both real and imaginary parts of A(q, cj) 
are important (Eq. 40). Because A'(q,cj) is even in cj, 
residual QPI intensity around the X symmetry point per¬ 
sists even at low energies, due to finite A'(q, cj = 0) = b'^. 
A'(q, cj) also diverges logarithmically (as Eq. 48) on ap¬ 
proaching the singular lines from X (it does not diverge 
in the vicinity of F or M). As such, the QPI scattering 
intensity is enhanced around X for magnetic impurities. 

Further, as shown in the lower panel of Fig. 3, the sign 
of Ap(q, uj) can be reversed by the contribution from the 
second term in Eq. 40. This is a hallmark of scattering 
from magnetic impurities, where Im Gd{oj) < 0 can be¬ 
come large due to the Kondo effect (see Eig. 2). This 
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leads to additional structure in the measurable |Ap(q)|, 
not found in QPI for scalar impurities. 


C. 3d cubic lattice with (100) surface 

The (100) surface of the 3d cubic lattice is again a 
square lattice; but ‘surface’ quasiparticles are dephased 
by coupling to the bulk. This leads to the t-matrix 
Eq. 24, and for a single impurity 

A(q,w)= / G°„rf(k||,w)G'°„rf(k|| - q,w) (49) 

ISBZ 

with scattering vectors q = qy confined to the 2d surface 
(and the Gg^^f(k||,cj) given by Eq. 25). 

The resulting QPI, shown in Eig. 4, does not contain 
the divergences arising in the 2d square lattice — but 
remnants of this singular structure appear in broadened 
regions of enhanced scattering intensity around the M 
symmetry point in the cubic lattice. The global four-fold 
symmetry of the QPI evolves locally into a continuous 
rotational symmetry around this point, with A(q, cj) = 
A(|q—qM|, cj) for |cj| <^t. Eurther, the QPI for the scalar 
impurity (which only depends on A"(q, cj)) is distinctly 
conical, with 

A"(q,w) aM(a;) +&m(w) |q - qivil + G |q - qivil^ , 

(50) 

whereas the QPI for a magnetic impurity can become 
dominated by the quadratic term (see Eig. 4) since 

A'(q,w) + cmH |q- qivil^ + G|q- qivil'^ • 

(51) 

A striking feature of the QPI for the 3d cubic lattices 
is the difference in intensity between scalar and magnetic 
impurities (note the rescaled color range in Eig. 4). There 
are two distinct reasons for this. Eirst, A"(q, uj) is odd in 
cj, whence A"(q, cj) ^ u; at low energies jcjl ^ t. Eor the 
scalar impurity, Eq. 38 implies that the QPI, Ap(q) ^ cj, 
is therefore also small. By contrast, the QPI for a mag¬ 
netic impurity (Eq. 40) has a contribution from A'(q, cj), 
which remains finite as cj ^ 0. 

Second, in the case of magnetic impurities, the Kondo 
effect produces a spectral resonance in Im Gd{oj) of width 
Tk that does not decay at low energies (cf. Eq. 43 and 
Eig. 2). In consequence, the QPI is considerably more in¬ 
tense at low energies for magnetic impurities than scalar 
impurities in standard flat-band metallic systems. 


D. 2d honeycomb lattice 

The 2d honeycomb lattice, modelling the tt and 
TT* bands in graphene within a nearest-neighbor tight- 
binding picture,generates richer structure in the 
QPI, 13,59,60 bipartite nature of the lattice and 

the low-energy pseudogapped density of states (Eq. 41). 




aoqcc/TT 



FIG. 5. QPI for a single impurity on the honeycomb lattice. 
As Fig. 3 but at uj = 0.3t and with impurity parameters from 
Fig. 2. The Brillouin zone cut takes the path F ^ M ^ 
F 2 ^ M 2 ^ K ^ F, where qr = 0, qM = |Ai, qps = Ai, 
qM 2 = Ai -h |A 2 , qK = |Ai -h |A 2 , in terms of reciprocal 
lattice vectors Ai ,2 = 27r/ao(^, ±1). 


A single impurity coupled to a single honeycomb site 
(on sublattice 7 ' = A or B) lowers the symmetry by 
breaking the centrosymmetry of the lattice. The single¬ 
impurity QPI is therefore obtained from Eq. 19 and 
Eq. 36, with A(q, cj) comprising contributions from both 
sublattices. 


A(q,a;) = ^Ajy^'(q,a;). (52) 

7 

A^’J (q, itself is given by Eq. 35, and depends on 
the phase 0 defined in Eq. 27b. This phase has a marked 
qualitative effect on the resulting QPI, and cannot be ne¬ 
glected. Eor centrosymmetric lattices, the QPI is periodic 
across the first Brillouin zone because A(q + nAi^uj) = 
A(q, cj) with integer n for any reciprocal lattice vector 
Ai. But in the honeycomb lattice 

A(q -I- SnAi.uj) = A(q, uj) , (53) 

arising because (/)(k + nAi) = 0(k) + exp[2n7ri/3].^^ As 
such, the period of the QPI is enlarged to include the 
third Brillouin zone. 

One consequence of this is the inequivalence of in¬ 
travalley scattering at the F and F 2 points (located at 
qr = 0 and qr 2 = Ai respectively), where Ap{qj^^,uj) = 
Ap(qr,ci;) x exp[—7ri/3]. At low energies |ci;| ^ t, these 
points are surrounded by singular lines in A(q, cj) at 
q = q*, where 


|q-q*| =dr(w) , (54) 

giving rise to a circular derivative discontinuity in the 
QPI around qr and qrs- The dispersive properties are 

















10 


controlled by dr{oo) = (ir 2 (<^), discussed further in the 
following subsection. For the scalar impurity at low en¬ 
ergies \uj\ ^ t, both points are surrounded by flat regions 
of high scattering intensity, 

l^p(q>w)l , (55) 

with br = br 2 independent of scanning energy oo. How¬ 
ever, the local environment of the F and r 2 points is dif¬ 
ferent. The immediate vicinity of the r 2 point possesses 
a continuous rotational symmetry, with divergences in 
the QPI along the entire singular line q = q*, 

A(q, cj) In^ |aoq — aoq* I . (56) 

By contrast, a lower six-fold symmetry is found around 
the F point as cj —> 0 due to divergent points arising only 
when q* X ^ = 0, with 8 = Ai, A 2 and (Ai + A 2 ). 

Complex features in the QPI also appear in the vicin¬ 
ity of the K symmetry points due to intervalley scat¬ 
tering, and are again enclosed by singular lines, de¬ 
noted q*. These features possess only reflection sym¬ 
metry about the line (qr — qx), the continuous rota¬ 
tional symmetry being lifted by the underlying phase 
texture (itself arising because the impurity couples to 
a single sublattice). The line of divergence along q* is 
intersected by a perpendicular nodal line at pinch-points 
where (q* — qx) * (qr — qx) = 0. For the scalar impurity, 
scattering is forbidden within the region around the K 
point enclosed by the singular lines. These features are 
seen clearly in the QPI map and cuts for the scalar im¬ 
purity presented in Fig. 5. At higher scanning energies, 
trigonal warping sets in, giving rise to a local three-fold 
point symmetry around qx- 

For the magnetic impurity, the relative weight of 
A'(q, cj) and A"(q, cj) in the QPI depends on the com¬ 
plex t-matrix, T(co’), which evolves with scanning energy 
uj (see Fig. 2). Importantly, this can lead to distinctive 
features in the measurable QPI, \Ap{q_,uj)\. Accidental 
cancellation of terms in Eq. 19 can produce ‘dark spots’ 
of suppressed scattering in the QPI. An example is shown 
Fig. 5, where |A/)(q, cj)! 0 for |q| < |q*| in the vicinity 
of the F point. In contrast to the scalar impurity case 
(Eq. 55), the QPI in general depends on uj and varies 
with q in the vicinity of the F and r 2 points when mag¬ 
netic impurities are present. Indeed, magnetic impurities 
also induce scattering near the K point. 


For scalar impurities, the scanning-energy dependence 
of the QPI is due entirely to the cj-dependence of A(q, cj), 
which characterizes the clean host lattice. The real part 
of this function is plotted as a colour map in the center 
column panels of Fig. 6 (the real and imaginary parts are 
related by Hilbert transformation). 

For magnetic impurities, QPI dynamics result from 
both A(q,cj) and the impurity Green function Gd{oj), 
whose spectrum is plotted in the left column panels of 
Fig. 6 (see also Fig. 2 and Eqs. 42-44 for the detailed 
low-energy behavior). The nontrivial scanning-energy 
dependence of the QPI reflects the rich structure of the 
underlying quantum impurity problem. 

For the 2d square lattice, divergences in A(q, cj), de¬ 
scribed by Eqs. 46 and 47, give rise to lines of intense 
scattering in the QPI. The dispersive properties of these 
features are controlled at low energies \uj\ ^ t by ar(co’) 
and aM(<^), which are related by continuity at the edge 
of the Brillouin zone through aui^) = ttr(co’)^/(27r). We 
find linear dispersion of the divergent features around the 
M symmetry point, implying 

aM(w) ~ |w| ; ar(w) ~ , (57) 

as confirmed directly in the upper panels of Eig. 6. 

As cj ^ 0, the divergences are confined to the line 
Qx ~ Qy connecting F and M symmetry points, whence 

■^(^Qx — Qyi ^ (2q^(C(j) 5q|C(j| T iCq T ... . (58) 

UJ 

The QPI itself thus diverges along this line with the uni¬ 
versal asymptotic form. 


\^p{qx = qy,(^)\ 


|u;|^0 


1 

(jj 


|u;| \n^ {IQt/\uj\) 


scalar 
mag. , 


(59) 


where Re Gd{j^ ^ 0) ^ [ln(16t/|co’|)]“^ is used in the case 
of the magnetic impurity (obtained by Hilbert transform 
of Eq. 43). The divergence is thus sharper along the M-F 
line for magnetic impurities, as evident in Fig. 6. 

Away from this divergent line (e.g. along the cut F ^ 
X ^ M in Fig. 6), the QPI is characterized by vanishing 
scattering intensity at low scanning energies due to 


HQxy^Qy,^) dq + 5q|cj|+iCqCjln|cj/t| + ... , (60) 

giving rise to the asymptotic behavior of the QPI, 


V. CHARACTERISTIC KONDO PHYSICS 
IN THE QPI 


|Ap(fe 7^ %,w)| 


|u;|^0 


Co’ln|co’/t| : scalar 
RWH7 • “""S- 


(61) 


A. Scanning-energy dependence 

We turn now to dynamical features of the QPI for 
the three lattices, comparing scalar and magnetic impu¬ 
rities. Numerically-exact results which exemplify the key 
physics are presented in Fig. 6 


As a result, the QPI intensity for the magnetic impurity 
decays much more slowly than that of the scalar impurity 
as the scanning energy is reduced in the vicinity of the 
X symmetry point — see Fig. 6. 

At higher scanning energies \uj\ ^ I//2, the ‘Hubbard 
satellites’ in the spectral function of the magnetic impu¬ 
rity give rise to enhanced scattering near the M and X 
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FIG. 6. QPI dynamics for a single impurity on various lattices: 2d square, 3d cubic with (100) surface, honeycomb (top, middle, 
and bottom row panels respectively). Right column panels: colour plots of the QPI |Ap(q,ct;)| across a Brillouin zone cut as 
a function of scanning energy cj/t, comparing scalar and magnetic impurities. Center column panels show a colour plot of the 
host function A'(q, u) over the same q-cut and energies. Left column panels show the spectral function for a magnetic impurity, 
—lm[tGd{^)\ vs cj/t, calculated via NRG at T = 0. Magnetic impurity parameters as in Fig. 2. Scalar impurity v = 0.15t. 


points. These features are not of course present in the 
scalar impurity QPI, and are as such one signature of 
strong electron correlations in magnetic impurities. 

For the 3d cubic lattice the dynamics are rather differ¬ 
ent, for two reasons: the host function A(q, cj) does not 
contain divergences, and the magnetic impurity Green 
function does not vanish at low energies because the host 
density of states is essentially flat for \uj\ <^t. 

For a scalar impurity on the cubic (100) surface, the 
QPI intensity vanishes everywhere at low energies |ci;| ^ 
0 because A"(q, uj) = — A"(q, —uj) is odd in u; — see cen¬ 
ter row panels of Fig. 6. By contrast, the Kondo effect 
gives rise to enhanced scattering at temperatures and en¬ 
ergies <C Tx] this gives rise to a large finite QPI intensity 
for UJ Tk- This is the typical behavior expected for 
magnetic impurities in standard metallic systems. 

Finally, we consider QPI dynamics on the the hon¬ 
eycomb lattice. In the case of the scalar impurity, the 
region of intense intravalley scattering around F and F 2 
described by Eq. 55 disperses linearly at low energies ac¬ 
cording to Eq. 54, with dr{uj) ^ uj. As cj ^ 0, the only 
divergent points are at q* = qp and qr2 • At q = qx , the 
QPI Ap(qK , cj) = 0 is identically zero for any |u;| < t. 
These features are shown in the lower panels of Eig. 6. 

The QPI for the magnetic impurity near these points 


shows intense scattering at energies |ci;| ^ U/2, corre¬ 
sponding to the Hubbard satellites in the impurity spec¬ 
tral function due to charge fluctuations. At lower ener¬ 
gies, however, Gd{oj) vanishes linearly according to Eq. 44 
due to the host LDOS which also vanishes linearly at low 
energies. Importantly, the Kondo effect is suppressed at 
particle-hole symmetry, and the local moment phase is 
always stable for any interaction strength.Electron 
correlations give rise to the nontrivial spin-flip scattering 
typical of such degenerate non-Eermi liquid phases. 

As a result the QPI for the magnetic impurity vanishes 
everywhere at low energies, according to 


|A/3(q,a;)| 


| u ;|^0 


|ci;| In^ \^/t\ 

< 0 

|ci;| In \uj/t\ 


q = qr(2) 

q = qx (62) 
elsewhere . 


B. Universality 

We focus now on a single magnetic impurity on the 
(100) surface of a 3d cubic lattice — the case most rele¬ 
vant to standard metallic systems where the host den¬ 
sity of states becomes essentially flat at low energies 
(see Eq. 41). The Kondo effect is operative in such 
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systems,^ with a spectral resonance setting in on temper- 
ature/energy scales ^ Tk — see center panel of Fig. 2. 
This resonance embodies enhanced spin-flip scattering, 
which screens the impurity local moment dynamically. 
Importantly, all physical properties depend only on the 
single emergent scale Tk at low temperatures/energies, 
reflecting the universal RG flow between local moment 
and strong coupling fixed points.^ 

At sufficiently low temperatures T ^ T^, the impu¬ 
rity Green function Gd{oj) is a universal function oiujjTK 
on all energy scales |c(;| ^ min(t, — not only for 

|cj| ^ Tk where strong coupling Fermi liquid behavior 
Eq. 42 holds, but also for \uj\ ^ Tk where local moment 
physics dominates the dynamics. In that case, the impu¬ 
rity spectral function takes the asymptotic form,^^ 


-Im Gd{uj) 




1 

1 + aln^ \buj/TK\ ’ 


(63) 


with a^b = 0(1) constants. This behavior for Tk ^ 
jcjl ^ min(t, is universal because the hybridiza¬ 

tion function Im T{uj) is essentially constant for |cc;| ^ 
min(t, V‘^/U) on the cubic lattice (and Re r(co’) ^ uj). 

Similarly, the real part of the host function A'(q, uj) be¬ 
comes constant on energy scales |cc;| ^ t, while A'(q, uj ) ^ 
UJ. In the scaling limit Tk 0, Im Gd{oj) thus controls 
the energy-dependence of the QPI in the universal regime 
(see Eq. 40). In practice, non-universal contributions are 
negligible for finite Tk ^ min(t, V‘^/U). 

In consequence, the entire QPI develops a universal 
scanning-energy dependence at low energies and tem¬ 
peratures. This means that magnetic impurities with 
different interaction strengths U and couplings V on 
different metallic substrates give the same normalized 
low-temperature/energy QPI | Ap(q, ci;)|/| Ap(q, 0)| when 
plotted vs rescaled scanning-energy for any scat¬ 

tering vector q. Rescaled experimental QPI data for dif¬ 
ferent systems should thus collapse onto a part of this 
universal curve, providing the unambiguous signature of 
scattering from magnetic impurities. This scaling col¬ 
lapse is demonstrated for the 3d cubic lattice in Eig. 7. 
Departure from universality is governed by the onset of 
ci;-dependence in A(q, cj), which is distinct for each q; 
for the experimentally relevant parameters used, this is 
found to occur at |ci;| > 20Tk — O.It. 


C. Thermal effects 

So far we have considered T = 0, appropriate in prac¬ 
tice when T ^ Tk ^ t, such that temperature is the 
smallest energy scale in the problem. This regime is 
relevant, as typical STM experiments are conducted at 
^ 5K;^^ but studies at higher T may also be performed^^ 
to investigate the change in QPI upon increasing temper¬ 
ature through Tk, and beyond. 

Non-interacting conduction electrons and uncorrelated 
impurities (e.g. the scalar impurity) in practice have T- 
independent electronic structure, and thus QPI. By con- 



FIG. 7. QPI I Ap(q, ct;)|/| Ap(q, 0)1 vs uj/Tk for a magnetic 
impurity on the 3d cubic (100) surface for scattering vectors 
Q = Qr, qx and qM (vertically offset by 0.5 for clarity), plot¬ 
ted for a range of impurity parameters: V/t = 0.555, 0.474, 
0.392 with fixed U = I.95t, corresponding to Tx/t = 5 x 10“^, 
7 X 10“^ and 3 x 10“^. Eq. 42 (solid line) and Eq. 63 (dashed) 
describe |ct;| Tk and ^ Tk asymptotes. 



aoqa;/7r Q^oq.x/T^ acq^A 


FIG. 8. Simulated experimental FT-STS measurements 
IApmeas(q,Cc; = 0 .ITk)| for a single magnetic impurity (pa¬ 
rameters as in Fig. 2, such that Tk = 5 x I0“^t) at fixed 
scanning energy, plotted across the IBZ for a series of tem¬ 
peratures, T/Tk = O.I, 0.2, 0.5, I, 2 and 5. 


trust, electronic correlations of the magnetic adatom ex¬ 
hibit a strong T-dependence, entering the t-matrix via 
the impurity Green function. As T is increased in metal¬ 
lic systems, the Kondo singlet is broken, destroying the 
Kondo resonance on the scale T > Tk and resulting in 
local-moment physics.^ This results in a quite dramatic 
change in the t-matrix, and hence QPI, on increasing T 
through Tk- 

In addition to this interaction-driven T-dependence, 
the local tunneling current measured in STS is weakly T- 
dependent due to thermal excitation of conduction elec¬ 
trons; at flnite-T, the differential conductance 

dl 

— {ri,Lo = eV,T)(x J ds p{ri,e,T)f{Lo - e,T), 

(64) 

where f'{uj - e,T) = ^/(w - e,T) and f{x,T) = [1 + 
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FIG. 9. QPI IAy9meas(qM)| for a single magnetic impurity 
(parameters given in Fig. 2, such that Tk = 5 x at the 

point qm = (tt, tt) in the IBZ, shown as a function of a; (T) 
at a series of different T (cj) in the upper (lower) panel; such 
that T/Tk (uj/Tk) = 0.01,0.1,0.2,0.5,1,2,5,10. Increasing 
T from low to high is denoted by the arrow. The QPI due to 
a scalar impurity is also plotted (dashed line) for comparison. 
Lower panel inset: the measured FT-STS (red) and the ‘pure’ 
QPI (blue), vs T/Tk, showing the same universal behavior in 
each case (up to a scale factor due to thermal broadening). 

exp(x/T)]“^ is the Fermi function. Eq. 64 represents 
the convolution of the LDOS (T-dependent only for the 
magnetic impurity) with a broadening kernel, controlled 
by T. The QPI measured via FT-STS, 
is then related to the ‘pure’ QPI by: 

/ oo 

ds Ap{q,s,T)f{u) - e,T) (65) 

-OO 

Fig. 8 shows the thermal evolution of the QPI for an 
impurity embedded on the 3d cubic (100) surface. For 
the magnetic impurity, the magnitude of the QPI de¬ 
creases substantially as T increases through Tk, and the 
Kondo resonance is suppressed. By contrast the QPI for 
a scalar impurity has a much simpler T-dependence (en¬ 
tering only via the thermal broadening, Eq. 65), with es¬ 
sentially no T-dependence for T ^ t, as depicted by the 
dashed line in Eig. 9. The strong T-dependence of QPI 
is a characteristic signature of Kondo physics in systems 
with magnetic impurities. 

The asymptotic low-T expansion of the impurity Green 
function at particle-hole symmetry and cj = 0 (applicable 
in the regime uj <^T ^ Tk) is a universal function of the 
Kondo temperature,^ 

Imr(w = 0) X Im Gd{uj = 0, T) 1 - ariT/T^^if + ... 

( 66 ) 

This leads to a close correspondence between the uj- and 
T-dependence of the QPI (via the impurity Green func¬ 
tion, comparing Eqs. 42,66), as seen by comparison of 
the upper and lower panels in Eig. 9. 


As such we expect to observe universal scaling in the T- 
dependence of the QPI for the magnetic impurity, analo¬ 
gous to that of the cj-dependence. The comparison of the 
pure QPI and the thermally broadened ET-STS signal in 
the inset of Eig. 9 (inset) demonstrates this universal be¬ 
havior, which is unaffected by the thermal broadening of 
the STM-measured conductance Eq. 65 (up to a trivial 
scale factor). 

VI. INTERPRETATION OF FT-STS 
A. Finite-size effects 

In experiment, the surface LDOS Ap{vi^uj) is mea¬ 
sured over an T X T plaquette using STM,^^ with the 
QPI obtained from Eq. 15. The q-space resolution of the 
resulting QPI naturally depends on the real-space sam¬ 
ple size. The ‘true’ QPI is recovered as T ^ oo, obtained 
theoretically by the t-matrix approach. 

We now consider explicitly the effects of finite sample 
size, by simulating the experimental protocol. The LDOS 
for these surface sites is calculated exactly using Eq. 16, 
with the non-local free Green functions obtained using 
the convolution method described in the appendix. 

Eig. 10 shows a Brillouin zone cut through the QPI 
for a magnetic impurity on the cubic lattice (100) sur¬ 
face, computed in the L ^ oo limit using the t-matrix 
approach (solid line). This true QPI is compared with re¬ 
sults for the same system restricting to an T x T surface 
sample, with L = 100 (crosses) and T = 10 (diamonds). 
The true QPI is very well-approximated when L = 100 is 
used (corresponding to a plaquette of side length ^10^ 
A, as typical in experiments^). The q-space resolution is 
also sufficient to capture accurately all features. Indeed, 
even for an extremely small sample region T = 10, the 
accuracy is surprisingly good; although the discretization 
is severe. 

Reassuringly, the experimental protocol reproduces ac¬ 
curately the true QPI. A large sample size is however still 
needed to resolve sharp q-space features; in 2d systems, 
the characteristic sharp cusps in the QPI in Eigs 3 and 5 
would require very large LDOS samples in real space. 

B. JDOS interpretation 

The QPI Ap(q, uj) is often interpreted (e.g. Refs. 40,13) 
in terms of the joint density of states (JDOS, J^(q, cj)), 
viz Ap(q, cj) J^(q, cj) with 

/)°(k, w)/9°(k - q, w) (67) 

ISBZ 

and p^(k,co’) = —G^(k,cj) the k-resolved surface 
density of states at energy uj. Contributions to the JDOS 
arise when quasiparticles on a constant energy contour 
with momenta k and k' are separated by q = k — k'. 
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FIG. 10. Brillouin zone cut of | Ay9(q, ct;)| for a magnetic impu¬ 
rity on the 3d cubic (100) surface at a; = 10“^t. Exact QPI 
(line) calculated via the t-matrix approach (Eq. 40), com¬ 
pared with the direct real-space approach (Eq. 15), sampling 
an L X L surface plaquette with L = 100 (cross points) and 
L = 10 (diamond points). Impurity parameters as in Eig. 2. 


On heuristic grounds, it is usually argued that the am¬ 
plitude of impurity-induced scattering from k to k' at 
energy uj can only be significant if there is a high density 
of quasiparticle states at both k and k' (i.e. they have a 
large JDOS). The QPI at scattering vector q is the sum 
of all scattering processes where k' = k — q. By assuming 
that the QPI is large when the JDOS is large, QPI pat¬ 
terns can be used to infer the JDOS and hence electronic 
properties of the clean host material. Experimental QPI 
data are typically interpreted in this way.^^’^^ 

The usefulness and relative simplicity of the JDOS pic¬ 
ture has motivated efforts to connect rigorously the QPI 
and JDOS. Notably, the perturbative approach employed 
in Ref. 64 makes the link by assuming a constant scatter¬ 
ing amplitude and phase along the constant energy con¬ 
tour. But a faithful description of interfering scattering 
processes typically requires relative phase information; 
and the JDOS simply lacks information about overlap 
matrix elements between states in the impurity-coupled 
system. In consequence, the QPI may be small even when 
the JDOS is large (as may be verified explicitly). The 
JDOS picture then fails to capture the basic physics of 
the scattering — as is known e.g. in graphene.We em¬ 
phasize that J^(q, cj) cannot be derived from Ap(q, cj) 
at a given q and cj, in any controlled limit. 

However for a simple scalar impurity on a centrosym- 
metric lattice, the QPI and JDOS are in fact related by 
Hilbert transformation. From Eqs. 38 and 45 (or e.g. 
Eq. 49), the QPI in this case can be written as 

Ap(ci,uj)=2v j [Re G°(k,w)]/)°(k-q,a;) , 

ISBZ 

( 68 ) 


where we have exploited periodicity across the ISBZ. We 


Ap(q.^' 


27rt;J'^(q, uj) 




the JDOS (right column) plotted as 27ru J^(q, a;) (see Eq. 69); 
for the 2d square lattice (upper row) and (100) surface of 
3d cubic lattice (lower row), for scanning energy u and Born 
impurity scattering potential v as in Eig. 3 and 4 respectively. 


introduce the complex quantity J(q, cj), defined as 
1 C d^k 

>^(q,‘^) =-/ — G'°(k,w)/)°(k-q,w) , (69a) 

TT J \IbZ 

ISBZ 

= - ;AAp(q, w) + i(q, w) , (69b) 

IV'K 

such that the QPI and JDOS are Hilbert conjugates, be¬ 
ing respectively the real and imaginary parts of J(q, cj). 

The JDOS interpretation of the QPI may thus be 
roughly correct for dilute scalar impurities on centrosym- 
metric lattices, because the underlying q-space topology 
of singular and nodal lines is the same for J^(q, cj) as it 
is for Ap(q, cj), both being related to a single complex 
analytic function J(q, cj). Nevertheless, even in this case 
the QPI is not accessible directly from the JDOS at a 
given q and uo: the entire cj-dependence of J^(q, uj) must 
be known to obtain Ap(q, uj) by Hilbert transformation. 
By way of illustration. Fig. 11 compares the JDOS to the 
QPI for the 2d square and (100) surface of the 3d cubic 
lattice; in the former case the two quantities are roughly 
similar, while in the latter the JDOS has significantly 
different q-space structure and intensity. 

We add that the above connection (Eq. 69b) does not 
hold in the case of non-centrosymmetric lattices or bi¬ 
partite lattices, where Eq. 68 is inapplicable. The QPI, 
Ap(q, cj), then depends on both real and imaginary parts 
of J(q, cj) due to additional q-dependent phase factors. 
This explains the failure of the JDOS interpretation in 
the case of impurities in graphene^^ (see Sec. IV D). 

Finally, we emphasize that the JDOS and QPI are 
never related simply by Hilbert transformation when 
magnetic impurities are present. This is because the 
t-matrix is a complex dynamical object: the real and 
imaginary parts of J(q, cj) are again mixed. Indeed the 
Hilbert transform, involving integration over all energy 
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scales, necessitates a full knowledge of the impurity dy- a program of the Netherlands Organisation for Scientific 
namics — information simply not contained in the JDOS. Research (AKM). 


VII. CONCLUSION 


APPENDIX: CONVOLUTION METHOD FOR 
LATTICE GREEN FUNCTIONS 


We have studied theoretically the use of quasiparticle 
interference (QPI) — measured in FT-STS experiments 
— as a probe of magnetic adatoms on surfaces. Fol¬ 
lowing a general formulation of the QPI due to an arbi¬ 
trary distribution of impurities, we turned explicitly to 
single-impurities adsorbed on a range of host surfaces: 
the (100) surface of a 3d simple cubic lattice, and the 2d 
honeycomb and square lattices, in which the Fermi-level 
densities of states respectively embody standard metallic 
behavior, pseudogap behavior, and a divergence due to 
a Van Hove singularity. In all cases, the single-impurity 
QPI factorizes into a local scattering t-matrix, and a host 
response function A(q, cj) at scanning energy uj = eV. 

The scattering t-matrix for a magnetic impurity is sim¬ 
ply related to the impurity Green function — itself de¬ 
pendent on the host lattice — and thus the rich dynamics 
due to electronic correlations is manifest in the QPI. De¬ 
spite the local, momentum-independent nature of these 
correlations, the q-space structure of the QPI is found 
to be qualitatively different from that of a simple scalar 
impurity due to non-trivial phase shifts associated with 
scattering from magnetic impurities, which reflect e.g. 
the Kondo effect. 

The response function A(q, cj) is also non-trivial, de¬ 
spite being a property of the free, non-interacting host. 
It displays significant structure in q-space, symptomatic 
of the symmetry and dimensionality of the host, but 
its energy-dependence becomes featureless for jcjl ^ t 
(with t the intersite lattice hopping). By contrast, the 
Kondo physics due to a magnetic impurity is controlled 
by an emergent scale Tk ^ t, and so scattering becomes 
strongly energy-dependent at low energies. Indeed, the 
QPI exhibits universal scaling in terms of cj/Tk and 
T/Tk — a characteristic hallmark for systems contain¬ 
ing magnetic impurities. Conversely, the QPI for systems 
containing scalar impurities has no energy or tempera¬ 
ture dependence on scales ^ t. 

The more complex case of QPI for multiple, mutually- 
interacting magnetic impurities, remains to be investi¬ 
gated. Such systems will display an even wider array 
of impurity physics, due to the competition of local and 
non-local (RKKY-type) interactions between impurities. 
These are expected to have a significant impact on the 
QPI for randomly distributed impurities, and will be the 
subject of future work.^^ 
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The calculation of real space lattice Green functions 
(LGEs) for periodic tight-binding (TB) models is a well- 
known problem relevant to many areas of physics.^^’^^’^^ 
The (retarded) real-space Green function between sites 
at r and r' on a d-dimensional lattice is generally given 
by, 


G°{r,r',oj)= I 

IBZ 


d^k 

^BZ 


^-i(r'-r).k 


Lu -|- i0“^ ~ ^0 ~ 


(70) 


where is the dispersion and cq = /i is a constant onsite 
energy (or chemical potential). Eor a nearest neighbour 
(NN) TB lattice specified by Eq. e^ = t with 

{dn} the set of NN lattice vectors. 

In Id, simple expressions for the local (on-site) and 
non-local (inter-site) LGEs are readily obtained, either 
by direct evaluation of Eq. 70, or via equations of 
motion.Eor the terminal site of a semi-infinite Id 
chain with on-site energies cq, the local LGE is given ex¬ 
actly in closed form by. 


G'?d(eo, w) = / where 


tf{uj) =U}- 


sgn(a))\/a)2 — 1 

iVl — 


|w| > 1 
|w| < 1 


(71) 


which is equivalent to Eq. 25. The LGEs for the infinite 
Id chain can be obtained in terms of Eq. 71 by exploiting 
translational invariance. 


G%^{x,x' ,eo,uj) 


w - eo - 2eG'l^{eo,uj) 


(72) 


However, LGEs for various lattice geometries 
in two and three dimensions are typically highly 
complicated,and not available in closed form. 
Direct numerical evaluation of Eq. 70 is notoriously diffi¬ 
cult, particularly for low energies, large site separations, 
or in the vacinity of Van Hove singularities. Recursion 
relations have been established in several cases, but 
solutions are often numerically unstable.Improved 
variants of the recursion technique (or continued fraction 
expansions) have been developed,but are costly to 
implement if LGEs are needed as an entire function of 
frequency. 

Here we derive a novel approach to the calculation of 
LGEs on hypercubic-type lattices, which is both highly 
accurate and numerically efficient. The method exploits 
the simple closed form expressions for the LGEs in Id, 
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Eqs. 71 and 72, building up lattices in higher dimen¬ 
sions by successive convolutions of those functions. The 
process is highly efficient because fast Fourier transform 
algorithms can be used to perform the convolution inte¬ 
grals. The method also has the advantage that boundary 
edges in 2d systems or explicit surfaces in 3d systems can 
be simply treated. Indeed, 2d nanoribbon or 3d block 
geometries can be implemented; and infinite systems can 
also be handled directly with no extra cost. 

We demonstrate the method first for the infinite 2d 
square lattice (and henceforth set cq = 0 without loss of 
generality). We denote the creation operator at lattice 
site r = {x^y) as cj ^ = Cy l^ and define the vector of 
operators for row ^ as = (...,cJJ., c^J., c^J., ...). The 
TB Hamiltonian then takes the form. 



FIG. 12. (a) 2d square lattice 



FIG. 12. (b) 3d cubic lattice with (100) surface 


^ ^ ^ ^ \^y ,(j ^ ^^7/,cr^y+l,cr H.C.^ , 

y= — oo a 

(73) 

where Tid a matrix describing the connectivity between 
sites of the ^-th row (here equal to t for nearest-neighbor 
sites and 0 otherwise). Importantly, for hypercubic-type 
lattices, Ti^ is independent of row index y. Eq. 73 repre¬ 
sents a set of coupled infinite Id TB chains to form the 
2d lattice. 

We now perform a canonical transformation of opera¬ 
tors = U^Cy^^^ with the matrix U defined such that 

D = WTidU is diagonal {D^k' = ^k^^kk')- Since this 
system is infinite and periodic in the x-direction, k = kx 
can be understood as the Bloch momentum, and ejf the 
Id dispersion. However, in general (e.g. for systems with 
a boundary), k merely labels an eigenstate of [/, with 
eigenvalue 

In this basis, Eq. 73 reduces to 


'H2d = E 

k 


E ^ {fy,afy+l,a + H.c.) 


(74) 


which describes a set of decoupled Id chains labelled by /c, 
each with constant on-site energy e]f. The transforma¬ 
tion from coupled to uncoupled chains is shown schemat¬ 
ically in Fig. I2a. 

The LGFs can then be expressed as. 


G'2d(r,r',w) = ((c“ V 


k,k' 




(75) 


where the last line follows because fy, oc 6kk' 

is diagonal in k (see Fig. I2a). The Id Green function 
given by Eq. 72. 


FIG. 12. Schematic showing the transformation from the real- 
space basis (left) to a basis of decoupled Id chains (right). 
Bonds denote hoppings t connecting sites. IJ diagonalizes 
rows of constant Vy in the infinite 2d square lattice system 
of (a); while M diagonalizes planes of constant > 0 in the 
semi-infinite 3d cubic lattice system of (b). 


We now make use of the spectral representation of the 
Id LFGs, which can be expressed in terms of and Uxk^ 
viz: 

Im (x, x', 0, Lo') = -TT y] U^^Ul^,5{uj' - el^) . 

k 

(76) 


One can then write Eq. 75 as, 

1 

C2d(r,r',w) = -- doj' Im G°^^{x,x' 

^ J —OO 

xC?d^(y,y',w',w) . 

(77) 

Since G^^^{y,y',u)',u:) = G^^^{y - y',u) - w') from 
Eqs. 71 and 72, Eq. 77 takes the form of a convolution 
integral. Gonvolution theorem then allows the efficiency 
of fast Fourier transform algorithms to be exploited, as 

J’u;[G°(r,r',w)] = -lF^\\.mG\^^{x,x',Lo)] 

^:F^[Gl^Jy,y',w)], 

(78) 


where denotes Fourier transformation. 

Straightforward extension of this method allows ac¬ 
cess to LGFs in higher dimensions. As a final instructive 
example, we consider now the 3d cubic lattice with an 
explicit (100) surface. The Hamiltonian is written as, 

OO 

Um = E E - t {cz\aCz+l,cr + H-C.) , 

z=0 a 

( 79 ) 
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in terms of vectors of operators for planes stacked in the 
^-direction, (f^^J, d^/J, ...), where 

creates an electron at site r = (x^y) of plane z. T 2 d 
is now the connectivity matrix for the 2d square lattice 
planes. 


The surface LGFs, with z = 0, then follow as 


G°urf(r,r',w) ={{cly,CoJ))^ 


We now diagonalize each 2d plane by writing f^a — 
such that D = M^T 2 dM is diagonal. As before, 
Dkk' = ^k^^kk'^ but now eff is the 2d square lattice dis¬ 
persion. In the transformed basis, the semi-infinite 3d 
cubic lattice becomes a bundle of decoupled semi-infinite 
Id chains, each with on-site energy as depicted in 
Fig. 12b 


where is given by Eq. 71. Employing the 

spectral representation of the 2d square lattice Green 
functions, 

Im r', w) = -tt y] M^^M^^,S{u} - (82) 

k 

we can write 

Gsurf(r,r',w) =-/ dco' Im G^rf(r,r',w') 

^ (83) 



y-sd = E 

k,a 


E 1 (/bvhi,.+ H.c.) 


.z=0 


(80) 


which can again be viewed as a convolution, here be¬ 
tween the semi-infinite Id Green function = 

G^^{uj — uj') given in Eq. 71, and the 2d square lattice 
Green function G 2 ^(r, r', cj) given in Eq. 77. 
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